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ABSTRACT 

Context. The presence of magnetic fields in the crust of neutron stars causes a non-spherically symmetric temperature distribution. 
The strong temperature dependence of the magnetic diffusivity and thermal conductivity, together with the heat generated by magnetic 
dissipation, couple the magnetic and thermal evolution of NSs, that cannot be formulated as separated one-dimensional problems. 
Aims. We study the mutual influence of thermal and magnetic evolution in a neutron star's crust in axial symmetry. Taking into 
account realistic microphysical inputs, we find the heat released by Joule effect consistent with the circulation of currents in the crust, 
and we incorporate its effects in 2-dimensional cooling calculations. 

Methods. We solve the induction equation numerically using a hybrid method (spectral in angles, but a finite-differences scheme in the 
radial direction), coupled to the thermal diffusion equation. To improve the boundary conditions, we also revisit the envelope stationary 
solutions updating the well known Tj, - Ts-relations to include the effect of 2~D heat transfer calculations and new microphysical 
inputs. 

Results. We present the first long term 2-dimensional simulations of the coupled magneto-thermal evolution of neutron stars. This 
substantially improves previous works in which a very crude approximation in at least one of the parts (thermal or magnetic diffusion) 
has been adopted. Our results show that the feedback between Joule heating and magnetic diffusion is strong, resulting in a faster 
dissipation of the stronger fields during the first lO"" - 10* years of a NS's life. As a consequence, all neutron stars born with fields 
larger than a critical value (> 5x 10"G) reach similar field strengths (a 2-3xlO'''G) at late times. Irrespectively of the initial magnetic 
field strength, after 10* years the temperature becomes so low that the magnetic diffusion timescale becomes longer than the typical 
ages of radio-pulsars, thus resulting in apparently no dissipation of the field in old NS. We also confirm the strong correlation between 
the magnetic field and the surface temperature of relatively young NSs discussed in preliminary works. The effective temperature of 
models with strong internal toroidal components are systematically higher than those of models with purely poloidal fields, due to the 
additional energy reservoir stored in the toroidal field that is gradually released as the field dissipates. 
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1. Introduction 

■ The neutron star (NS) magnetic field (MF) maintained by elec- 
' trie currents circulating in the crust modifies the crustal temper- 
' ature distribution by means of two mechanisms. The first one 
is due to the anisotropy of thermal conductivity in presence of 
a strong MF and causes important changes on how heat flux 
flows from the core through the crust and envelope up to the 
surface. The second mechanism is the generation of heat due 
to MF dissipation which results in a non spherically symmetric 
allocation of heating sources. Therefore, the magnetic diff'usiv- 
ity, besides its tensorial character due to the presence of the MF 
(which eventually results in the Hall drift of the field), becomes 
a non-spherically symmetric quantity. Under these conditions, 
an initially purely dipolar field evolves very quickly in a magne- 
tized NS to generate complex structures including toroidal fields 
and/or higher order multipoles. Different studies on magneto- 
hydrostatic equilibrium configurations indicate that stable con- 
figu rations require both toroidal and poloidal components (see 
e.g. iReiseneggej (l2008h and references therein). 

The complexity of the problem has limited previous works 
to partial studies of the complete problem. The anisotropic heat 
flux and its consequences for the - in principle observable - 
surface tempe rature (Ts) distribution has been considered by 
iGeppert et all (|2004 12006) and iPerez-Azorfn et all (l2006allbb . 
but for fixed, prescribed MF. The long-term effect of Joule heat- 



in g has also been recently di scussed by lAguilera et al.l (I2008d) 
or lUrpin & Konenkovl (l2008h (see also references the r ein), al- 
though in a very crude approach. In lAguilera et al.l (I2008bl) 
a full 2-D cooling code has been used to describe the tem- 
perature evolution, but the Joule heating rate is estimated by 
an anal ytical approximati on and considered to be uniform. In 
Urp in & Konenkovl (l2008b the former results of iMiralles et al.l 
(1998|) are revisited in the context of high field radio-pulsars, 
but with the simplified approach of considering the diffusion 
of a one-mode (dipolar) poloidal field, without actually per- 
forming consistent simulations including the temperature evo- 
lution. On the other hand, the evolution of the MF including 
th e influence of the Hall term at early times has been studied 
m IFons & Gepperll (l2007l) . The Hall drift causes a somewhat 
faster dissipation for strong initial fields due to the reorganiza- 
tion of the field in smaller scales. In this latter work the temporal 
evolution of the temperature has been prescribed according to 
a generally accepted cooling law and assumed to be the same 
in the whole crust. The effect of ambipolar diffusion in the liq- 
uid core can also b e relevant, as recently studied for example in 
iHovos et al.1 (l2008l) . 

One of the reasons that leads to complex MF geometries in 
the crust is due to the non-spherically symmetric magnetic diffu- 
sivity. It is effective even if the Hall drift is negligible. Since the 
temperature within the crust is no longer uniform, the magnetic 
diffusivity in this approach becomes a function of the radial and 
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of the polar coordinate, thus rendering the one-mode approxi- 
mation inappropriate. Up to now, the effect of the dependence 
of the magnetic diffusivity on the polar angle (if axial symme- 
try is assumed) has not been considered. In some studies about 
cooling of NSs, the source term in the thermal diffusion equa- 
tion includes an angular-averaged Joule heating rate, where the 
heat production is uniform in spherical shells at a given radius 
r (Pageetal. 2000; Aguilera et al. 2008c; Urpin & Konenkov 
These are not fully consistent calculations, because the 
currents can locally be very intense and release heat in a highly 
non-spherically symmetric way. Therefore, it is important to at- 
tack the problem of the full magneto-thermal evolution of NSs, 
by solving simultaneously the induction equation and the heat 
transfer equation taking into account the anisotropy of the ther- 
mal conductivity and the electrical resistivity tensors in a consis- 
tent way. This is the main goal of this paper, in which we solve 
this problem (in axial symmetry) for the first time. The aim of 
this study is to provide a self-consistent model for the coupled 
evolution of the MF and temperature in axially symmetric NSs, 
in general relativity, and including the effect of angular varia- 
tions of temperature in the magnetic diffusivity and the thermal 
conductivity. We also revisit the envelope stationary solutions 
updating the well known T\, - Tj-relations to include the effect of 
2-D heat transfer calculations and new advances in microphysi- 
cal ingredients (phonons, ion-ion interactions). We provide new 
fits to solutions of equilibrium magnetic envelope that can be 
used as boundary conditions in multidimensional cooling simu- 
lations. 

The paper is organized as follows. In the next Section 
we derive and present the basic equations together with the 
boundary conditions. Next we describe the input microphysics 
and discuss the effect of superfluid neutrons on the heat flux 
anisotropy. The following Section is devoted to the presentation 
of the magneto-thermal evolution. Finally we discuss the results 
by comparing them with observable consequences. 



2. Basic Equations 

The thermal evol ution of the NS is c alculated using the 
code described in lAguilera et al.l (l2008bl) . For the evolution 
of the MF we b a sically follow the formalism already used in 
iPons & GeppertI (|2007|) . We address the interested reader to 
these references for more details about the cooling code and the 
MF evolution code that we have merged and coupled to perform 
the simulations presented later in this paper. In the next subsec- 
tions we just summarize the basic equations that are solved with 
the purpose of introducing notation and for the sake of complete- 
ness. 



2.1. Heat Transfer Equation 

Assuming that deformations with respect to the spherically sym- 
metric case due to rotation, MF, and temperature distribution do 
not affect the metric in the interior of a NS, we use the standard 
static metric 



-c'e 



(1) 



where v(r) is the relativistic redshift correction and A(r) is the 
length correction factor, e _ ^ _ 2GM{r)/c^r. Using this 
background metric, the thermal evolution of a NS can be de- 
scribed by the energy balance equation 



dT , 
ot 



where Cv is the specific heat per unit volume and Qy is the energy 
loss by neutrino emission while Qi, stands for the Joule heat- 
ing rate. In this study we neglect other heating mechanisms than 
Joule heating, such as accretion or frictional heating at the core 
crust interface. In the diffusion limit, the heat flux is simply 



-e-"k ■ V(eT) 



(3) 



where k is the thermal conductivity tensor that depends on both 
coordinates {r,6) through its dependence on temperature and 
MF. Note that the differential operator V associated to the spatial 
metric e^'^^''^dr^ + P'dQ} must include the corresponding metric 
factors (i.e. e"''^''^^ for the radial component). An explicit rep- 
resentation in terms of the metric scale factors can be found in 
Radler et al. (2001). We solve Eqs. ^ a nd Q numeric a lly fol - 
lowing the same scheme a s described in iGeppert et al] (|2004|) ; 
iPerez-Azonn et all ( l2006al) . 

2.2. The Induction Equation 

The Maxwell equations (in the Gaussian system) relative to an 
observer at rest (Eulerian observer) for the static metric given by 
Eq. ([T]i are 

IdE „ , V „^ 47r 
_ = Vx(e''B)- — 

c ot c 

VB = 

= -Vx(e''£) (4) 

c at 

where again the differential operators are defined according to 
the spatial coordinate system, E and B are the electric and mag- 
netic field measured by the Eulerian observer, as well as and 
j that represent the electric charge and current density, respec- 
tively. 

In the comoving frame, which corresponds to the Eulerian 
frame since matter is at rest, electric field and density current are 
related through Ohm's law. This law establishes a proportional- 
ity between the current and the electric field which can be written 
as j - crE, cr being the electric conductivity. In the presence of 
strong MFs, cr becomes a tensor and Ohm's law must be written 
as j - &E or equivalently E - 'Rj, where = o" ' is the resis- 
tivity tensor. This te nsor can be decomposed i n symmetric and 
antisymmetric parts (" Landau & Lifshitall960l) . The symmetric 
part, for non-quantizing fields is isotropic and determined by the 
scalar — , where crw is the electric conductivity in the direction of 

o-|| II 

the MF. The antisymmetric part of the resistivity tensor can be 
represented by a vector proportional to B. The electric field E is 
then written in the form (IZimanll 19791) 



E = —j + -^B X j 

cry erigC 



(5) 



where is the electron number density and e is the elementary 
charge. 

If we neglect the displacement current term in the Ampere- 
Maxwell equation (second equation in |4]i, and make use of Eq. 
(|5]l, the Faraday induction equation (fourth equation inUJi can be 
written as 

^ = -V X X (e^'B) + [V X (e^'B)] x b\ (6) 



(2) where we have introduced the magnetic diffusivity t] = 
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There ar e two main differ e nces f rom the induction equation 
employed in iPons & Gepperj (l2007h : the relativistic factors are 
included and the magnetic diffusivity is not assumed to be spher- 
ically symmetric. 

Since both tie and the metric factor only depend on the ra- 
dial coordinate (it is well justified to neglect the structural defor- 
mations induced by MFs) the non-line ar Hall term can be treated 
in the same way as done by iPons & G eppert (2007). However, 
the magnetic diffusivity through its temperature dependence also 
depends on the polar angle, — T](r, 0), thus requiring an exten- 
sion of the previous formalism. Th e influence of th e Hall term 
at early times has been studied in P ons & GeppertI (l2007h . re- 
sulting in a somewhat faster dissipation for strong initial fields 
due to the reorganization of the field in smaller scales. The nu- 
merical treatment of the non-linear term is complex and very 
computationally-limited. In this paper our goal is to perform 
coupled (magneto-thermal) long term simulations, thus from 
now on we focus on the linear part of the induction equation. 

The hnear part of the induction equation (|6]l reads 



and 



^ = -V X [77 V X (e^'B)] . 



(7) 



We deco mpose the MF into its poloidal and toroidal part 
dRadleret al. 2001) 



B - Bpoi + Bt( 



(8) 



and use their representation in terms of two scalar functions 
<l)(r, e) and "fir, 0): 



B 



pol 



-V X (r X V (D) , fit. 



-rx V^P. 



(9) 



We expand now the functions <1>, ^, and ?; in a series of spherical 
harmonics (in the axisymmetric case) as follows: 



1 °° 
r ^ 

n=\ 

= -V »P„(r,Oi'„(0) 

00 

77 = Yj^„{r,t)YM, 



(10) 



where Y„ is the spherical harmonic K'" for m - 0. 

The poloidal and toroidal parts of the MF can be written as 



B 



pol 



= — ^ n(n + 1)<1>„F„ e,- + — ^ 



" d0 



dY„ 
dr d0 



(11) 



where e^, eg, are unitary vectors associated to the spatial met- 
ric. Inserting this expressions into Eq.|2]and using the quantities 
7^,^ and defined in Geppert & Wiebicke (1991, Eqs. 59 and 
60) which are related to the Clebsch-Gordan coeflicients, we 
arrive at the following evolution equations for the poloidal and 
toroidal scalar functions: 

k.k' 



-2A 



dv dX\ d^k k(k + 1) 



[dr drj dr 



(12) 



k(k + 1) 



k'k 



(13) 



If 77 = ri(r, t) only k' - contributes and, taking the non- 
relativistic limit (e" - = 1), the purely diffusive Eqs. 8 of 
iPons & Geppert! (I2OO7I) with D„„, - C„,n - are recovered. Note 
that no coupling between poloidal and toroidal component ap- 
pears, since this coupling can only be obtained by including the 
Hall drift term. 



2.3. Microphysics 

The microphysical ingredients that enter in the heat transport 
and induction equation are the specific heat, the thermal con- 
ductivity, the neutrino emissivity and the magnetic diffusivity. In 
the solid crust, the dominant contribution to the specific heat is 
that from electrons and ions, while electrons, lattice phonons and 
collective modes of superfluid neutrons contribute to the ther- 
mal conductivity. We refer the reader to a detailed description 
of the employed equations of state, the thermal conductivity, the 
spe cific heat, and the ne utrino emissivities, given in section 4 
of iAguileraetal.1 ilOOm. In this pap er we assume the mini- 
mal cooling scenario dPage et aLll2004l) which includes neutrino 
emission from the Cooper pair breaking and formation process, 
but we do not take into account any direct Urea process, due ei- 
ther to nucleons or to exotic matter (hyperons, Bose condensates 
or deconfined quarks). 

In addition, it has been shown by IChugunov & Haensell 
dlOOTt) that the ionic contribution to the total thermal conduc- 
tivity is negligible in the crust but it may play a role for low tem- 
peratures in the envelope. Very recently dAguilera et alj|2008al) 
studied the possible effects of collective modes of superfluid 
neutrons. They found that this process may dominate the thermal 
conductivity in the inner crust when its temperature is ^ lO^K. 
For such relatively old NSs, heat transport by superfluid neutrons 
counteracts the anisotropy in the electron conductivity caused by 
a strong crustal field and, eventually, turns the inner crust isother- 
mal. In this work we have updated our microphysics and we have 
also included this two new contributions to the thermal conduc- 
tivity. 

The only relevant contributions to the electrical resis- 
tivity are electron-p honon and electron-impurity collisions 



dFlowers & Itoh|[T976l) . While the efficiency of electron-phonon 
collisions strongly depends on the crustal temperature, the 
electron-impurity scattering is much less sensitive to it. 
However, both processes may be strongly affected by the pres- 
ence of a strong MF which suppresses the conductivity compo- 
nents perp endicu lar to the field lines ( |C anuto & C hiuderi 1970| ; 
Iltohlll975l) . As in|P ons & Geppert! d2007l) . we calculate the elec- 
trical resistivity by using the electron relaxation time provided 
by A. Potekhin's public codeQ We have used an impurity con- 
centration parameter of 0. 1 (the definition of this parameter and a 
discussion about how it affects electronic tran s port can be fo und 
in section 5.1.1 of lPerez-Azorm et al.l (l2006a!) ) . iJonesI dl999h has 
shown that disorder in the inner crust could result in an impurity 
parameter > 10, which leads to larger electrical resistivity and 
enhanced ohmic decay. However, this effect becomes important 
only for sufficiently cool NSs. 



www. iof f e . rssi . ru/astro/conduct/condmag . html 
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2.4. Joule Heating 

Joule heating couples the thermal and magnetic evolution by 
contributing to the source term Q/, in Eq. |2] which, for large 
fields, can result in a higher temperature and therefore a larger 
magnetic diffusivity. These feedback may lead to a faster dissi- 
pation of the MF if it is strong enough to really alter the tempera- 
ture of the crust. In this paper Joule heating is taken consistently 
into account in a cooling simulation for the first time. As the MF 
evolves in time, we compute at each point of the computational 
grid the local values of the electrical current density, which is 
simply 



j(r,0,t) = e""— Vx(e''B) . 
4-n 

The corresponding components of the current density are 

47T 1 v-l 

c 

n 



(14) 



z 



-X 



-2A 



dr de 
d^^n 2GM(r) (90,, n(n+l) 



d0 



We evaluate the heating source term by 

,■2 



(15) 



(16) 



where Qh is the energy per unit time and unit volume measured 
by the Eulerian observer. The magnetic energy balance equation 
is 



d 



dtV Stt 



'—ExB 

4-n 



(17) 



which is easily interpreted: if we integrate over the whole vol- 
ume, the magnetic energy losses are due to Joule heating and 
Poynting flux through the boundaries. This latter term vanishes 
since we do not consider the possibility of having electromag- 
netic waves (we neglect displacement currents). 

2.5. Magnetic boundary conditions 

Since we restrict ourself to MF configurations confined to the 
crust, the inner boundary conditions are determined by the re- 
quirement that the normal component of the MF and the tangen- 
tial components of the electric field has to vanish at r = Re. This 
is a consequence of the assumption that the core is in a super- 
conducting state and the Meissner-Ochsenfeld eff'ect prevents 
the MF to penetrate. Therefore we apply the following boundary 
conditions atr-Rr 



dr 







,2 



k.k' 



(18) 
(19) 



A detailed derivation is given in iPons & Gepperj (l2007b . Note 
that in the limit of vanishing Hall-drift this reduces to 



dr 



= 0. 



For the outer BC we require all components of the MF to 
be continuous across r - 7?ns to match the relativistic vacuum 
solution. Hence, let us first consider the stationary solution for 



the outer space. In the absence of external currents, the toroidal 
component of the MF must vanish, and each multipole of the 
poloidal field must satisfy (in the stationary case); 



(1 -z)-r^ + -- 



n(n + 1) 



dr 



(20) 



where z - IGMjc^r, which corresponds to the compactness pa- 
rameter at r = /?NS- 

This second order diff'erential equation has analytical solu- 
tions for each value of n, although they cannot be written in a 
closed analytic form valid for any n. For example, for n - \ and 
n = 2 we have 



Oi = Cir" 



O2 = Cir^ 



ln(l -z)+z+ — 



Cir^ 



(4-3z)ln(l -z)+Az-z^-- 





(21) 



+ C2r3(4-3z)(22) 



where C\ and C2 are arbitrary integration constants that must be 
fixed according to the value of the magnetic multipole moments. 
Regularity of the external solution at r = 00 requires C2 = 0. 
In general, the family of solutions of Eq. |20] for any value of 
n can also be expressed in terms of generalized hypergeometric 
functions {F{{], [],z)), also known as Barnes's extended hyper- 
geometric functions, as follows: 



O,, = Cir-" F({n,n + 2],{'^ + 2n\,z) 

+ C2r"+'F([l-«,-l-n],[-2n],z). 



(23) 



Note that regularity at r = 00 requires again C2 = 0. 

For practical reasons, we write the outer boundary conditions 
at r = 7?NS as 



dr 



fn^n 







(24) 
(25) 



where /„ is a relativistic factor that only depends on the com- 
pactness ratio z{Rm) (in the Newtonian limit /„ = 1) and can be 
evaluated numerically or with the hel p of any algebraic m anipu- 
lator using the form given in Eq. dia. Radleretal.l(l2001h use an 
alternative form based on the expansion of the vacuum solutions 
in a series of powers of 1/r. 

3. Thermal boundary conditions: blanketing 
envelopes revisited 

The very diff'erent thermal relaxation timescales of the envelope 
and the crust of NSs makes computationally expensive any at- 
tempt to perform cooling simulations in a numerical grid that in- 
cludes both regions simultaneously. Since radiative equilibrium 
is established in the low density region much faster than the crust 
evolves, the usual approach is to use results of stationary, plane- 
parallel, envelope models to obtain a phenomenological fit that 
relates the temperature at the bottom of the envelope Tb, with 
the surface temperature T^. This Ts - Ts{Tb) phenomenologi- 
cal function is used to implement boundary conditions, because 
it allows to calculate the surface flux for a given temperature at 
the base of the envelope. is generally chosen to correspond 
to some density between the neutron drip point p ^ 3 x 10"g 
cm"-' and p - lO'^g cm"-'. Examples of such models of magne- 
tized e nvelopes have been c onstructed by [P otekhin & YakovlevI 
(1200 ll) and later upgraded in lPotekhin et al.l (i2007.) to include the 
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Fig. 1. Surface temperature as a function of temperature at the 
neutron drip point for two models of NS envelopes with different 
MF strengths. The thin dash ed lines are the analytical fits from 
jPotekhin & Yakovlevll200ll) for MF parallel (||) and perpendic- 
ular (±) to the normal direction to the surface. The thick dashes 
are our results for the ± case, and the solid lines are our full 2D 
results for a dipolar field (we show at the pole and at the equa- 
tor). Our results for a purely radial field cannot be distinguished 
from the upper solid line (J ^ at the pole for the dipolar MF). 
Numerical results always include the effect of neutrino emission 
in the outer crust and envelope. 



effect of the neutrino emissivity in the outer crust. They derived 
an analytic form of the Tb - relation that reads 



r,(B, ^, g, Tb) ^ Tf\g, Tb) V, Tb), 
where 



(26) 



(27) 



and ^ = 0.in,g - 0.001 g^^^ ^jQ.l Tb,&. Here gi4 is the surface 
gravity in units of 10'"* cm s"^, Tb.s is in 10** K, and r^.g is 
in lO*" K. 



The function X has been fitted by decomposing into transver- 
sal and longitudinal parts as 

XiB, ip, Tb) = [ X'^^l^iB, Tb) cos^ ,p + X^^iB, Tb) sin^ ipf^ ,(28) 

wher e (p is the angle between the MF and the normal to the sur- 
face. [Pote^in&^^tovle^lllOOl]) give the following fits for Xj^ 
and X\\ 



30.292j-0.240 



<Y||(B, Tb) = 1 -H 0M92B\f 

VI + 0.1076Bi2(0.03 + rfe.9)-"-559 



XAB, Tb) 



[1 +0.819Bi2/(0.03 + rfc,9)] 



0.6463 



(29) 
(30) 



where B12 = B in units of lO'^ G. These fits are valid for B < 
10'* G and 10^ K < Tb < K. 

It must be reminded that the above T\,-Ts relation is based on 
a plane-parallel approximation. When this approach is applied 
to a spherical star, meridional heat fluxes in the envelope are 
not allowed and, therefore, this approximation may be inaccurate 
when the se fluxes c ompete with the purely radial ones. In addi- 
tion, Chugurioy & Haensel (2007) found that the contribution of 
ions or phonons to the thermal conductivity of the envelope can 
reduce the anisotropy of heat conduction. Therefore, we have re- 
visited the magnetized envelope problem with two motivations: 
upgrading the microphysical inputs (thermal conductivity) and 
assessing on the accuracy of the plane-parallel approximation. 

To avoid solving the hydrostatic equilibrium equations in 
two dimensions, we have build a spherically symmetric iron en- 
velope model with a zero temperature equation of state. With this 
fixed background, we have calculated stationary solutions of the 
heat transport equation in 2D, with a given MF geometry (dipole 
solution). At first glance, this may seem inaccurate, since finite 
temperature effects may be relevant at low density but, as ex- 
plained in Gudmundsson et al. ( 1983), the main regulator of the 
Jb - Ts-relation is the sensitivity strip where the opacity is max- 
imum. This strip marks the transition from electronic heat trans- 
fer to the radiation dominated one. It lies at relatively high densi- 
ties except for very low temperatures. Hence, structural changes 
at low density do not affect th e solution, even for e xtreme cases 
such as a condensed surface jPotekhin et al.l [2007b . To test the 
validity of this assumption and to compare with previous works, 
we have taken the crustal MF to be either purely radial (labelled 
II, parallel to the normal to the surface) or purely tangential to 
the surface (labelled ±). In this two limiting cases, the heat flux 
is purely radial, a nd we have reproduced withi n a 5% the T\, - 
relation given bv lPotekhin & Yakovlevl(l2001h if neutrino emis- 
sivities are not considered. The minor differences are probably 
due to the different EOS or NS model (i.e. our neutron drip point 
is at p = 3.5 X 10" g/cm^, while theirs is at p = 4 x 10". Our 
results for purely parallel or purely tangential fields including 
the neutrino emissivity in the outer crust are show n with thick 
dashes in Fig. [T] The analytical fit from (Potekhin & YakovievI 
l2001h . without the effect of neutrino emissivity, is plotted with 
thin dashes. The agreement is very good in those cases where 
neutrino emission is not important (low temperature). Moreover, 
our results for radial and tangential MF are directly compara- 
ble to Figs. 4 and 5 of Potekhin et al. (2007), in which neutrino 
emissivity effects are considered. The very good agreement with 
this work indicates that the use of a zero temperature EOS for 
solving the hydrostatic equilibrium equation is a valid approxi- 
mation. The same comment applies to the effects of strong MF 
in the low density EOS, as long as they are only noticeable at 
densities where radiation dominates the heat transport. Finally, 
in order to check the effect of using an EOS in which the thermal 
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contribution to the pressure is neglected to obtain the mechanical 
structure of the envelope, we have also recalculated some of the 
models using new envelope structures obtained by solving the 
hydrostatic equilibrium equations for a finite temperature EOS 
and assuming the temperature profile of the previous models. 
We have found the same surface temperatures in both models 
(within a 1 %), further confirming the validity of our approach. 
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Fig. 2. Stationary 2D solutions of magnetized envelopes: surface 
temperature profiles as a function of the polar angle for Bp - 
10'^ G and three different combinations of the r^-distribution 
and MF geometry. The solid lines correspond to the numerical 
results and the dashed lines to the analytic fits from Eqs. (|3TT i 
and(l28]l. 

Concerning 2D models, solid lines in Fig. [1] correspond to 
our 2D transport results with a dipolar field. The quoted value of 
the MF corresponds to the strength at the magnetic pole. From 
our results, we confirm that the effect of ion/phonon transport in 



the envelope is to reduce the large anisotropy obtained in pre- 
vious magnetized models. In addition, by performing 2D heat 
transport calculations through the magnetized envelope, we take 
into account the meridional heat fluxes driven by the meridional 
temperature gradients between pole and equator, and we find that 
the anisotropy is further reduced. This effect is more relevant for 
high fields, resulting in an equatorial temperature about a fac- 
tor 3 lower than the polar temperature, in contrast with the 2-3 
orders of magnitude obtained in previous models. For practical 
purposes, we have made fits of our results keeping a similar func- 
tional form to Eq. [30l but only changing the form of X\\ and 
as follows: 

<Y||(B, Tb) = 1 + 0.05B°f 



VI + 0.07Bi2(0.03 + Tb,9)-"-^^'^ 
[1 + 0.9Bi2/(0.03 + Tbs)] 



0.4 



(31) 



We find that this formula reproduces very well the correct nu- 
merical results for the T\, - relation for a dipolar field geom- 
etry. The ratio between polar and equatorial temperature is sig- 
nificantly lower than in the case of 1-D plane-parallel envelope 
models. We also confirm that, if neutrino emission at densities 
p > lO'" g/cm^ is taken into account, there is a MF-dependent 
upper limit to the effective surface temperature (iPotekhin et al.l 
I2007 '). but our maximum temperatures for MF tangential to the 
surface are considerably higher than those obtained before. We 
find that the upper limit consistent with our numerical results is 
well fitted by the following expressions 



T^max 



3.6 X lO^'K 
1 -I- 0.02 log B 12 
2.8 X 10*'K 



1 +0.61ogBi2 



(32) 



It must be stressed that these upper limits are obtained by 
the inclusion of the outer crust (3 x 10" > p > 10'" g/cm-') 
in the static envelope calculations that fix the internal bound- 
ary condition. If the numerical grid of the cooling code extends 
to densities lower than 10'" g/cnr' this effect is consistently in- 
corporated and there is no need to include these corrections in 
the T\, - Ts relation. Note that we neglect the possibility of 
Joule heating within the envelope because, due to the very large 
magnetic diffusivity, initial electric currents are quickly dissi- 
pated and may only be important at the very beginning of a 
NSs life. Obviously these results may vary for other field ge- 
ometries. We have checked, however, that our improved fits ( [3T1 ) 
combined with Eq. ( l28T l are a good approximation (within 5%) 
to the numerical result in a number of cases. In Fig. ^ we 
compare the numerical results with the analytical fits for three 
different models. In the top panel the MF geometry is purely 
dipolar and the temperature distribution at the inner boundary 
is fixed to Tt = ^^(cos^fl -i- 0.2sin2 6') K. The middle panel 
shows the results for a model with the same MF configuration 
but a temperature distribution at the base of the envelope given 
by Th - 10^(1 + 3 sin^ 6) K. The bottom panel corresponds to 
a model with the same Ti, distribution as the middle panel, but 
with a purely quadrupolar MF geometry. In all cases Bp - 10'^ 
G, so that the combination of high field and low temperature 
demonstrates the validity of the fit for the most extreme cases. 

4. Results 

We present the results now of our numerical simulations of the 
magneto-thermal evolution of NSs and its dependence on the 
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20 kyr 



100 kyr 



500 kyr 





Fig. 4. Temperature distribution in the crust of a NS (i.e. up to the bottom of the envelope) at different ages. The initial MF is purely 
poloidal with {Bp = lO''* G). The field lines are also shown. The numbers in the color-scale on the left of each figure indicate the 
maximum and minimum values of the temperature (in units of 10** K) at each age. The crustal shell has been stretched a factor of 4 
for clarity. 




Fig. 3. Power spectrum at f = 5 x 10^ years for the models 
corresponding to purely poloidal fields with initial values of 
Bp = 10'^ (crosses), 10'^ (stars), and 10'^ G (diamonds). 



initial MF structure and strength. We restrict this pres entation 
to a sel ection of the initial models described in Ag uilera et al.l 
(I2008hl) . Our model A at f = consists of a dipolar (n = 1 ) 
poloidal field according to Eq. (13) of .Aguilera et al.l (l2008bl) . 
parameterized by the value of the radial component at the mag- 
netic pole. Bp. In model B we add to the initial dipolar field a 
quadrupo lar ( » = 2) toroidal component obeying Eq. (11) of 
[Aguilera et al.1 (l2008bl) . with a maximum value of x 50Bp. This 
latter model represents the case of a strong toroidal component 
confined in the crust. We recall that we focus in this paper on the 
purely diffusive case, i.e., the Hall term in the induction equation 
is neglected. This term needs a separate specific numerical treat- 
ment that is out of the current capabilities of our code and will 
be discussed in later work. Its influence can be important either 
at very early times or, as we point out at the end of the section, 
at late times during the photon cooling era. For most of the first 
10'' - 10^ years of a NS life we do not expect qualitative changes 
from our present results. 

To begin our discussion, we present in Fig.[3]the power spec- 
tra for model A with three initial fields of Bp = 10'^, lO''*, and 
10'^ G at the age of f - 5 x 10^ years. We see how the cou- 
pling between different modes due to the angular dependence of 



T] fills the shorter wavelength modes (initially only the dipolar 
poloidal component « = 1 is present). Only odd multipoles are 
present because the initial model is symmetric with respect to the 
equatorial plane. At this age, the cascade has filled out all large 
wavenumber modes and is saturated following approximately an 
n"^ power law for the cases with large initial fields. This result 
shows that, even if the influence of the non-linear Hall term is 
negligible, the thermo-magnetic coupled evolution results in a 
complex field geometry that can not be described by a single 
mode. Although the dominant mode is still the dipole, the dis- 
tribution of part of the energy in higher other modes leads to a 
slightly faster dissipation. 

In Fig.|4]we show the crustal temperature distribution of a 
NS at diff'erent ages (t - 10, 100 and 500 kyr), together with the 
poloidal field lines. The tendency of the surfaces of constant tem- 
perature to be aligned with t he magnetic field lin es discussed in 
iPerez-Azorin et al.l ( l2006al) : lGeppert et al] (12004) is clearly vis- 
ible, but there is an important difference with respect to previ- 
ous works. The consistent inclusion of the Joule heating source 
results in a greater energy deposition in the region where cur- 
rents are more intense. In this particular model this happens at 
low latitudes, since the currents that maintain a poloidal field are 
toroidal. Thus more heat is released close to the equator. While 
the polar regions are always in thermal equilibrium with the core, 
the equator is insulated due to the reduced thermal conductivity 
across magnetic field lines. As a consequence, the heat released 
at early times in the equatorial region cannot flow across field 
lines into the core, where it would be rapidly lost by neutrino 
emission. It is only allowed to flow towards the surface along 
field lines. This modifies the traditionally accepted temperature 
distribution consisting of hot poles and a cooler equator (insu- 
lated from the warmer core). Instead we find that the equatorial 
region of the outer crust is actually warmer than the poles. In 
models with weak MFs, when the effect of Joule heating be- 
comes less effective, the situation is inverted and the thermal 
distribution with hot polar caps is recovered. Notice the large 
difference between polar and equatorial temperatures, up to a 
factor of 2 at early times (the numbers next to the color scale 
indicate the maximum and minimum temperature at each evolu- 
tionary time), in contrast with the nearly isothermal crusts at late 
times. 

In non-magnetized NSs the heat conduction between the 
crust and the core is able to counteract the local neutrino emis- 
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sion processes and produce an isothermal NS after the thermal 
relaxation stage, which lasts 10-100 years. This interchange of 
energy between core and crust is still important in magnetized 
neutron stars, but the presence of a strong tangential field in the 
equatorial region of the inner crust has the effect of increasing 
the thermal relaxation time in several orders of magnitude, so 
that the equatorial region in Fig. |4] becomes effectively decou- 
pled from the core. 
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Fig. 5. Temperature profiles at the base of the envelope {Tt) and 
at the surface {T^) at different evolutionary times. The initial MF 
is purely poloidal with Bp - 10'"^ G. 

We must remind that the anisotropy in (at the bottom of 
the envelope), is not necessarily the same as that of T^. The blan- 
keting effect of the envelope, and/or atmospheric effects should 
be taken into account before a comparison with observations. 
In Fig. (|5]l we plot the evolution of both, the temperature at the 
bottom of the envelope (top panel) and the surface temperature 
(bottom panel) as a function of the polar angle. 

In Fig. |6l we show a sample of cooling curves (effective 
temperature versus true age) for models A (solid lines) and B 
(dashed lines) but varying the initial strength of the field. For 
high fields the effect of Joule heating is visible from the very be- 
ginning of the evolution. The effective temperature of a young. 




Fig. 6. Cooling curves. Effective temperature as a function 
of age for different initial field strengths (from bottom to top 
Bp = 10^^ 3 X 10'^ 10'4, 3 X 10l^ and 10'^ G). The solid lines 
correspond to Model A and the dashed lines to model B. 



t - 10^ yr magnetar with Bp - 10'^ G is higher by a factor 
of 2 than that of a NS with a standard Bp - lO'"* G, and it 
is kept nearly constant for a much longer time. The effect is 
further enhanced in model B due to the extra energy stored in 
the toroidal field. In spite of some quantitative differences, our 
i mproved simulat i ons co nfirm the qualitative results described 
in lAguilera et al ] d2008ch . where a phenomenological MF de- 
cay law has been adopted. The most important difference is the 
somewhat h igher temperature reac hed in these models, when 
compared to lAguilera et al ] (I2008ch . The difference stems from 
the particular location of the heat dissipated by the Joule effect 
in the regions where currents are intense, as opposed to the ho- 
mogeneous distribution of energy in the formerly adopted phe- 
nomenological model. 

In Fig.|2]we show the value of the MF at the pole as a func- 
tion of time for models A (solid lines) and B (dashed lines) with 
different initial field strengths. Naively, we can divide the mod- 
els in two groups: strong initial field Z? > 5 x 10'^ G and weak 
initial field Z? < 5 x 10'^ G. The plot shows that models with 
strong initial fields are subject to a faster decay than those with 
weaker fields. Models with weak initial fields are subject to de- 
cay in a more or less similar manner The field typically decays 
in about a factor of two on a timescale of few 10^ years, and then 
remains nearly constant, due to the increase of the magnetic dif- 
fusion timescale as the star cools down. 

Models with initially large fields behave in a different way. 
Since the magnetic energy stored in the crust is now large 
enough to significantly affect the thermal evolution when it is 
steadily released by Joule heating, it results in a higher average 
temperature of the crust. But this process has a back-reaction: 
the higher temperatures imply larger resistivities and therefore 
faster decay. Interestingly, at f > lO*" yr, all models with large 
poloidal initial fields seem to converge to an asymptotic fidu- 
4 X 10'^ G, while all models with strong 



cial value of B, 

toroidal fields converge to a lower value of B, 



asymp 



asymp 



1.5 X 10' 



G, because of the higher temperatures reached in average during 
the evolution. 

Based on these results we predict that all sufficiently old NSs 
born with Bp above a critical value evolve towards similar field 
strengths, while those born with lower fields show a MF distri- 
bution at late times similar to the initial distribution but shifted 
in about a factor of 2 to lower values. This critical initial field is 
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approximately B„h ~ 5 x 10'^ G. It also delimits the minimum 
field strength that can actually influence the thermal evolution 
of the NS by MF diffiision. Once the MF is dissipated below 
that value, its influence on the later evolution is reduced and the 
subsequent evolution proceeds in a similar way in all cases. Our 
models predict that no old magnetar can be found, and that, at 
ages f > 5 X 10^' yr, the typical MF for all NSs born as magnetars 
must be similar (B ^ Basymp)- Up to know, observational data are 
not in contradiction with this fact. 

To conclude this section we turn now to discuss the evolution 
of the dimensionless magnetization parameter cobt, where cjb is 
the electron gyro-frequency and r is the electron relaxation time. 
The Hall induction equation (|6]l can also be written as 

dB I \ \ 

^ = -^'^ X — X (^''^) + "^BT [V X ie^B)] X CB 1 , (33) 
ot An \(T|| / 

where en is the unit vector in the direction of B. 



100.0 ^ 



m 




0.1 I . . . 

10^ 10" 10^ 10^ 

t (yr) 

Fig. 7. Evolution of the MF strength at the pole Bp during the 
first million years of a NS for several initial values. The solid 
lines correspond to Model A and the dashed lines to model B. 

In this form, the interpretation of ojbt is straightforward: 
when the magnetization parameter exceeds unity, the Hall drift 
term dominates. The effect of the Hall drift term at early times 
has been discussed in detail in (Pons & Geppert 2007). As a rule 
of thumb, when cjbt > 10^ the effect of the Hall drift signifi- 
cantly changes the evolution, while more moderate values lead 
to a somewhat faster dissipation due to reorganization of the field 
in smaller scales. Due to numerical limitations, in this work we 
have restricted ourselves to the purely diffusive case, but it is 
worth to look at the evolution of this parameter as shown in Fig. 
[8] Here we show radial profiles of the magnetization parame- 
ter in the crust for model A and three different initial values of 
Bp. The upper panel corresponds to f = 10^ yr and the lower 
panel to f = 10^ years. At early times, the scaling of llibt with 
Bp is visible and only in magnetars one must expect large values 
in the inner crust. However, the back-reaction of the field evo- 
lution on the temperature has some interesting implication. At 
late times, the temperature of low field NSs is lower than that 
of highly magnetized NS, so that the temperature dependence of 
the electron relaxation time overcomes the effect of the MF and 
it turns out that the former magnetars are less magnetized than 
NSs born with moderate fields. In addition, at ages of 10^ yrs 
the temperature is low enough to ensure that the magnetization 
parameter is very large for all models studied, specially in the 



inner crust (where the conductivity is larger). This opens more 
questions about the late evolution of NSs that is likely to be dom- 
inated by the non-linear Hall term in most NSs. In particu lar, the 
possibility of the Hafl instability (iRheinhardt & Geppertii2002l) . 
or implications on the evolution of the braking index are worth 
to be explored. 




11.0 11.2 11.4 11.6 
R (km) 

Fig. 8. Radial equatorial profiles of the magnetization parameter 
cjbt for model A and three different initial values of Bp: 10'^ 
G (solid lines), lO''* G (dashed lines), and 10'^ G (dotted lines). 
The upper panel corresponds to f = 10"^ yr and the lower panel 
to f = 10^ years. 



5. Conclusions 

We have performed consistent 2D simulations of the coupled 
magneto-thermal evolution of NSs for the first time, including 
realistic microphysical input and general relativistic corrections. 
By properly taking into account the interplay of the MF and 
temperature evolution in cooling simulations we have found that 
their mutual influence is important and affects the outcome when 
the initial field strength is of the order or larger than a critical 
value of Bciit = 5 x 10'^ G. It has effects on both sides: the tem- 
perature and the field strength. 
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The average effective temperature of a NS born with Bp < 
Bcrit is barely affected, while those born as magnetars are sub- 
ject to significant heating by the dissipation of currents in the 
crust. In addition, since heating is locally important in the re- 
gions where currents are more intense, the surface temperature 
distribution can be very different depending on the field geome- 
try. While the pole (or other regions in which the field lines are 
nearly radial through the crust) is in thermal equilibrium with 
the core and has its same temperature, regions in which the field 
lines are tangential remain essentially thermally isolated. This 
can produce cooler areas if no heating process is considered or, 
conversely, hotter regions if heating is important. This is the case 
of poloidal fields, in which currents are located near the equato- 
rial region that remains warmer than the rest of the crust during 
a long time. The effective temperature of models with strong in- 
ternal toroidal components are systematically higher than that of 
models with purely poloidal fields, due to the additional energy 
reservoir stored in the toroidal field that is gradually released as 
the field dissipates. 

In the models with stronger fields, as a result of the aver- 
age higher temperatures, the crustal electrical resistivity is en- 
hanced and magnetic diffusion proceeds faster during the first 
10^ - 10'' years of a NS's life. As a consequence, all NSs born 
with fields larger than a critical value (> 5 x lO'-'G) reach simi- 
lar field strengths 2 - 3 x lO'-'G) at late times, irrespectively 
of the initial strength. After 10'' years the temperature is so low 
that the magnetic diffusion timescale becomes longer than the 
typical ages of radio-pulsars, resulting in apparently no dissipa- 
tion of the field in old NSs. We confirm the strong correlation 
between the MF and the surface temperature of relatively young 
NSs discussed in preliminary works. Notice that if the MF of 
magnetars is caused by superconducting currents in the core (as 
opposed to crustal currents) the longer diffusion timescale in the 
core would allow magnetars to live much longer with their orig- 
inal large fields. Thus, observations of magnetars can help to 
discern between models with currents located in the crust or in 
the core. 

It should be mentioned that magnetic field evolution by am- 
bipolar diffusion in the core may produce qualitatively similar 
effects to those obtained in this work: keeping the temperature 
high while the field is stron g and stopping fiel d decay when the 
temperature drops (see e.g. lReiseneggeij|2008l) . In this work we 
have not considered the MF evolution in the core because am- 
bipolar diffusion is usually considered under the assumption that 
the core is in a non-superfluid state and therefore is important 
during the very early stages of evolution. We are more interested 
in the long-term evolution, after the temperature rapidly drops 
below the critical temperature for nucleon superfluidity. In or- 
der to quantify the relative importance of both effects one would 
need to consider ambipolar diffusion in a superconducting fluid 
coupled to the dissipation of crustal currents studied here. 

The detection of magnetars with true ages (the spin-down 
age can be seriously overestimated) t a \{f yr, or the detection 
of a young highly magnetized NS with T < 10^ K would be a 
serious challenge for crustal field models. At present our results 
are in agreement with the known population of high field NSs 
and magnetars and support the idea of the existence of a strong 
crustal MF component in magnetars. 

Given the complexity of the feed-back between temperature 
and MF, it seems necessary to extend this work in two main lines 
that can shed new light on our knowledge of the cooling theory 
of NS: engaging 3D simulations and including the Hall term in 
the induction equation. The complex geometry that may arise in 
a realistic case with hot spots and irregular fields is certainly not 



treatable with our present code and needs further investigation. 
Similarly, the influence of the Hall term at late times (when the 
temperature is so low that the magnetic diffusivity is negligible) 
or a consistent treatment of ambipolar diffusion in a supercon- 
ducting fluid coupled to the dissipation of crustal currents, can 
produce new interesting effects and are issues worth to explore 
in future works. 

Acknowledgements. We thank D.N. Aguilera for valuable comments and pro- 
viding updated conductivity routines used in the simulations, and Andreas 
Reisenegger for a critical reading and his constructive and helpful comments. 
This research has been supported by the Spanish MEC grant AYA 2007-67626- 
C03-02 and the Research Network Program Compstar funded by the ESF. U. 
Geppert thanks the University of Alicante for support under its visitors program. 



References 

Aguilera, D. N., Cirigliano, V., Pons, J. A., Reddy, S., & Sharma, R. 2008a, 

ArXiv e-prints, 807 
Aguilera, D. N., Pons, J. A., & Miralles, J. A. 2008b, A&A, 486, 255 
Aguilera, D. N., Pons, J. A., & Miralles, J. A. 2008c, ApJ, 673, L167 
Canuto, V. & Chiuderi, C. 1970, Phys. Rev D, 1, 2219 
Chugunov, A. I. & Haensel, P 2007, MNRAS, 381, 1143 
Flowers, E. & Itoh, N. 1976, ApJ, 206, 218 
Geppert, U., Kiiker, M., & Page, D. 2004, A&A, 426, 267 
Geppert, U., Kiiker, M., & Page, D. 2006, A&A, 457, 937 
Gudmundsson, E. H., Pethick, C. J., & Epstein, R. I. 1983, ApJ, 272, 286 
Hoyos, J., Reisenegger, A., & Valdivia, J. A. 2008, A&A, 487, 789 
Itoh, N. 1975, MNRAS, 173, IP 
Jones, P B. 1999, Physical Review Letters, 83, 3589 

Landau, L. & Lifshitz, E. 1960, Electrodynamics of Continuous Media (Addison 
Wesley) 

Miralles, J. A., Urpin, V., & Konenkov, D. 1998, ApJ, 503, 368 
Page, D., Geppert, U., & Zannias, T. 2000, A&A, 360, 1052 
Page, D., Lattimer, J. M., Prakash, M., & Steiner, A. W. 2004, ApJS, 155, 623 
Perez-Azon'n, J. E, Miralles, J. A., & Pons, J. A. 2006a, A&A, 451, 1009 
Perez-Azorin, J. E, Pons, J. A., Miralles, J. A., & Miniutti, G. 2006b, A&A, 459, 
175 

Pons, J. A. & Geppert, U. 2007, A&A, 470, 303 

Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 2007, Ap&SS, 308, 353 
Potekhin, A. Y. & Yakovlev, D. G. 2001, A&A, 374, 213 

Radler, K.-H., Euchs, H., Geppert, U., Rheinhaidt, M., & Zannias, T. 2001, 

Phys. Rev D, 64, 083008 
Reisenegger, A. 2008, ArXiv e-prints, 809 

Rheinhardt, M. & Geppert, U. 2002, Physical Review Letters, 88, 101103 
Urpin, V. & Konenkov, D. 2008, A&A, 483, 223 

Ziman, J. 1979, Principles of the Theory of Solids (Cambridge University Press) 



